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ABSTRACT 

Nonaxisymmetric, meridional circulation inside a neutron star, excited by a glitch 
and persisting throughout the post-glitch relaxation phase, emits gravitational radia- 
tion. Here, it is shown that the current quadrupole contributes more strongly to the 
gravitational wave signal than the mass quadrupole evaluated in previous work. We 
calculate the signal-to-noise ratio for a coherent search and conclude that a large glitch 
may be detectable by second-generation interferometers like the Laser Interferometer 
Gravitational- Wave Observatory. It is shown that the viscosity and compressibility of 
bulk nuclear matter, as well as the stratification length-scale and inclination angle of 
the star, can be inferred from a gravitational wave detection in principle. 
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1 INTRODUCTION 

Rotation-powered radio pulsars are promising sources of high frequency gravitational waves. Their spin frequencies often lie in 
the hectohertz 'sweet spot' of current detectors, e.g. the Laser Interferometer Gravitational- Wave Observatory (LIGO). The 
rotation of their crusts can be measured extremely precisely, enabling coherent searches which improve the signal-to-noise 
ratio by the square root of the number of wave cycles obser ved. Such coheren t searches have already beaten el ectromagnetic 
spin- down limits on the quadrupole moment of the Crab ([Abbott et al. 2008) and are close for other pulsars ( Abbott et al. 
2007). There are two main obstacles to detection. (1) Dephasing occurs if the radio pulses are used to construct a grav- 
itational wave phase model but the fluid interior rotates at a slightly different speed to the crust. (2) The quadrupoles 
predicted so far are relatively small in isolated pulsar s without any ongoing a ccretion activ ity, e.g. unstable oscillatio ns such 
r- modes ([Brink et al.1 2004; iNavvar fe Owen 2006; Bondarescu et al. 2007), precession (j Jones &; Ande rsson 2002), inter- 



nal magnetic defo rmations ( Bonazzola &; Gourgoulhonl 1 19961: ICutlei 



Siderv et al.ll200d ). and hydrodynamic turbulence (Melatos & Peralta 201 



20021) qu 
l l20ld). A 



quasiradial fluctuations ([Sedrakian et al.l I2OO3I ; 
Ac creting millisecond pulsars can reach larger 



quad rupoles through magnet ically confined mountains jMelatos fe Payne] 20051 : Payne &; Melatos! l200d : Ivigelius &; Melatos] 
20091 ) or thermal mountains (jUshomirskv et alJbood : Iflaskell et alJbood ). 

In this paper, we investigate another source of gravitational radiation from i solated pulsars, namely the radiation emitted 
during the recovery phase following a pulsar glitch (jvan Eysden &; Melatos 20081 ) . Glitches are small, abrupt jumps Ais in the 



rotation frequency v which range in fractional size from 10 
individu al objects. Current ly, out of ~ 1800 known pulsars 



11 to 10 4 across the pulsar population and over four decades in 
101 have been observed to glitch, with a total of 285 individual 
events ([Melatos et al.ll2008|). Glitches occur ra ndomly in all but two objects (PSR J0537-6910 and PSR J0835-4510), which 



spin up quasiperiodicallv ([Melatos et al. 20081 ) . Most pulsars which have glitched at all have only glitched once. Of the 35 per 



cent that have glitched multiple times, and with the exception of the quasiperiodic p air, the glitch sizes a nd waiting times 
are well fitted by power-l aw and Poissonian probability density functions respective ly ( Melatos et al. 20081 ) . consistent with 

an avalanche mechanism (Warszawski & Melatos! [iooi : iMelatos fe Warszawskill200d ). 

Most theories of pulsar glitches build on the vortex unpinning paradigm introduced by Anderson &; Itohl (1975). Super- 
fluid vortices pin to lattice sites or defects in the crust and are prevented from migrating outward as the crust spins down 
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electromagnet ically. At some stage, many vortices unpin catastrophically, transferring angular momentum to the crust. While 
it is unknown what triggers the collective unpinning, it is likely to excite a nonaxisymmetric flow for two generic reasons. 
(1) Pinning causes the crust and superfluid to rotate differen tially, inevitably driving nonaxisymmetric merid i onal circulation 
and e ven turbulence, as observed in laboratory experiments (Munson fe Menguturk 1975; Nakabavashi 1983; Junk fe Egbers 



2000) and numerical simulations (IPeralta et al.ll2005l . 12006 al ibi: Melatos fe Peraltall2007l : IPeralta et al.ll2008l : IPeralta fe Melatos 



20091 ) of spherical Couette flow. (2) Avalanche trigger mechanisms, like self-organized criticality, which are favoured by the 



observed glitch statistics, intrinsically lead to an inhomogeneous and hence nonaxisymmetric superfluid velo city field, with 
spatial fluctuations correlated on all scales, from the smallest to the largest ( Jensen! 19981 ; Melatos et al. 2008 ). 

The gravitational wave signal from a pulsar glitch separates into two parts. First, there is a burst corresponding to 
nonaxisymmetric vortex unpinning and rearrangement during the spin-up event itself. To date, observations have failed to 
resolve th e spin-up time-scale. In the Vela pulsar, w hich was monitored continuously for several years, it occurs over less 
than 40 s ([McCulloch et al.lll990l : bodson et alll2002h . Second, there is a decaying continuous-wave signal during the quasi- 
exponential relaxation phase (lasting days to weeks) following the spin- up event ( Shemar fe Lv ne 1996). The latter signal arises 
as viscous interactions between the crustal lattice and core superfluid erase the nonaxisy mmetry in the sup erfluid velocity 



field and restore the crust and core to co-rotation (or at least steady differential rotation). Side rv et al.l (|2009l ) constructed 



two-fluid 'body- averaged' model of a glitch and calculated that the burst signal emitted during the spin- up event by coupling 
to quasiradial oscillations is too weak to be detected. In this paper, we focus on the second part of the signal, which has the 
advantage of enduring for many rotation periods, enabling a coherent search with increased signal to noise. 

Two techniques h ave been proposed to date to search for gravitational radiation emitted during the spin-up event and 
post-g litch relaxation . IClark et al.l (|20Q7h developed a Bayesian selection criterion for comparing f-mode ringdown to white 
noise. Havama et al. ( 20081 ) investigated coherent network analysis, which does not assume any particular waveform. Both 
methods would be aided by the availability of a specific signal template, like the one calculated in this paper. Importantly, 
by combining such a template with data, gravitational wave experiments can constrain the equation of state of bulk nuclear 
matter, complementing particle accelerator experiments which have recently produced results that disagree with astrophysical 
data. Heavy ion and nuc l ear resonance experim ents measuring the compressibility of nuclear matter imply a soft equation of 
state ( Sturm et al. 200ll; IVretenar et all 2003), whereas neutro n star observations imply a hard equation of state, albeit at 



lower energies ( Hartnack et al. 2006; Latti mer fe Prakashl 2007 ). Likewise, heavy-ion colliders measure a viscosity close to the 



conjectured quantum lower bound (Adar e et al ]|2007^ w hereas the relaxation t i me-scale of pulsar glitches s uggests a value 
many orders of magnitude larger (Cutler fe Lindblom 1987; Andersson et al . 20051; Ivan Evsden fe M elatos 201ol ). Gravitational 
wave observations will help to resolve these and other issues; bulk matter at nuclear density cannot b e assembled in terrestrial 
laboratories with current technology ( van Evsden fe Melatos! 2008 ; Owen et al. 20091 ; Xu et al. 2009). 

In this paper, we calculate the gravitational radiation generated from the spin up of the stellar interior following a 
pulsar glitch. We estimate its detect ability with the current generation of long-baseline interferometers, and show that certain 
important constitutive prop erties of a neutron star can be extracted from gravitational wave data, at least in principle. 
The calculation is based on van Evsden fe Melatos ( 20081 ) , extended to treat current quadrupole radiation. In Section [2] we 
solve the general hydrodynamic problem of nonaxisymmetric, stratified, compressible spin-up flow in a cylinder, driven by 
Ekman pumping , following an abrupt in crease in the angular velocity of the container. The initial and boundary conditions 
implemented by van Evsden fe Melatos (2008) are modified slightly to make them more realistic. In Section [3] we predict 
the gravitational radiation emitted during the relaxation phase following a glitch. We calculate the signal-to-noise ratio and 
estimate the detectability of the signal in Section 2] In Section [5] we show how to extract the compressibility, stratification, 
and viscosity of the stellar interior from gravitational wave data. 



2 EKMAN FLOW FOLLOWING A GLITCH 



Radio pulse timing experiments h ave so far failed to resolve temporally the abrupt increase in the angular velocity of the 
neutron star crust during a glitch ( McCulloch et al. 1990l : Dodson et al. 2002). Hence, in the absence of more detailed infor- 



mation, we model a glitch as a step increase in the angular velocity Q, of a rotating, rigid, cylindrical container filled with a 
Newtonian fluid ( Abnev fe Epstein! 19961 ; 



Evsden fe Melatos! 20081 ). A cylinder is a coarse approximation to a spherical 



star, but it admits analytic solutions and has a long history of being used to model neutron stars and in geomechanical studies 
( Pedloskvl[l967l : IWalinl(l969l : lAbnev fe Epstein! 1 19961 ; Ivan Evsden fe Melatos! 1 2008^ . 



Differ ential rotation betwee n the container and interior fluid drives Ekman pumping, which spins up the interior over 
time; see lBenton fe Clark (1974) for a review of Ekman pumping. The spin up of an axisymmetric container was first treated 
analytically bv lGreenspan fe Howard! (|l963l ). For an incompressible fluid, the entire volume is spun up on the Ekman time- 
scale, tE — i^ -1//2 n -1 , where E — v/(QL 2 ) defines the dimensionless Ekman number in terms of the kinematic viscosity v 
and the size L of the container. Subsequently, it was shown that compressibility and stratification reduce the spun-up volume 
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by hindering flow along the side walls (|Walinlll969l : Ubnev fe Epstein! 1 19961 : Ivan Evsden"fe Melatos 2008h. With less volume 
to spin up, the Ekman time-scale is lower. Nonaxisymmetric spin up was analysed by van Evsden &; MelatosI ( 2008h . 

In this section, we solve the problem of the nonaxisymmetric, stratified, compressible spin up of a cylinder, extending 
van Evsden &; Melatos ( 20081 ). We write down the linearised hydrodynamic equations in Section 12.11 solve for the general 



spin-up flow in Section 12.21 apply initial and boundary conditions in Section 12.31 and 12.41 and discuss precisely how and why 
these conditions differ from previous analyses. The final, time-dependent solutions for the pressure, density and velocity fields 
are presented in Section [2. 51 We discuss the initial condi t ions f or a glitch in Section [2. 61 For full details of the calculation, the 
reader is referred to Section 2 of van Evsden fe MelatosI (|2008h . 



2.1 Model equations 

Consider a cylinder of height 2L and radius L, containing a compressible, Newtonian fluid with uniform kinematic viscosity z/, 
and rotating about the z axis with angular velocity ft — £le z . In the rotating frame, the compressible Navier-Stokes equation 
reads 

v+t).Vu + 2ftxu = --Vp + g + vV 2 v + (V • v) + V ( Itfr 2 ) . (1) 
ot p 3 V2 J 

The fluid satisfies the continuity equation 
d f t + V • (pv) = , (2) 
and the energy equation is written in a form that relates the convective derivatives of the pressure and density, 

(i + " v )'=?(i + *- v >- (3) 

The symbols v, p, p, g and c represent the fluid velocity , density, pressure, gravit ational acceleration, and the speed of sound, 
which is determined by the equation of state. Following I Abnev &; Epstein! (1993), gravity is taken to be uniform and directed 
towards the midplane of the cylinder, 

-ge z if z > , 

+ge z if z < . { } 

where g is constant. 

We work in cylindrical coordinates (r,(/>,z) and consider the region z ^ 0, as the flow is symmetric about z = 0. Equations 
{J)-© are rewritten in dimensionless form by making the substitutions t H> t^t, r H> Lr, z H> Lz, v H> LSQv, p \-> pop, 
p 1— »• pogLp, and V 1— > (1/L)V, where the sca le factor pp is chosen to be the equilibrium density at z = 0. The scaled equations 
obtained in this way [see equations (6)-(8) in Ivan Evsden &; MelatosI (2008)] feature three dimensionless quantities: the Rossby 
number e = £0/Q, the Froude number F — LO 2 /g, and the scaled compressibility K — gL/c 2 . 



2.2 Spin-up flow 



At time t = 0, the angular velocity of the cylinder accelerates instantaneously from O to O + SQ. If e is small, as in a pulsar 
glitch, the problem linearises and we can solve for the equilibrium and spin-up flows separately by making the perturbation 
expansions p \-> p + e5p, p \-> p + e5p, and v \-> Sv. In the frame rotating at O, the equilibrium velocity is zero and the spin-up 
flow is of order e. 

We assume the equilibrium state is steady and axisymmetric, with p = p(r, z) and p = p(r, z). Ignoring centrifugal terms 
proportional to F, and takin g p~ 1 dp/dz to be uniform for simplicity, as in previous work l|Walinll 19691 : 1 Abnev fc Epst ein 1996; 
van Evsden fc MelatosI 12008). we find 



p{z) 



K7 1 e~ K " 



(5) 
(6) 



where K s = L/z s = —Lp~ 1 dp/dz is a constant which depends on the stratification length-scale, z s . 

The spin-up flow is uns teady and nonaxisymmetric, w ith Sp = Sp(r, 0, z, £), 5p — Sp(r, 0, z, t), and Sv = Sv(r, 0, z, t). We 
solve equations (17)-(21) in Ivan Evsden &; MelatosI (|2008l ) for the spin- up flow using the metho d of multiple scales, expanding 
Sv, 6p and dp as perturbation serie s in the small parameter E 1 /2 , e.g. Sp 5p° +E 1/2 5p 1 +Q {E) (jWalin 1969: 1 Abnev fe Epsteinl 
19961 : Ivan Evsden & MelatosI 1 20081 ). Following Section 2.3 in Ivan Evsden fc Melatoi <|2Q08h . the O(E ) continuity equation is 
automatically satisfied and the order 0(E 1 ^ 2 ) equations reduce to 



ld_ 

r dr 



(•• 



dr) 



+ 



1 d 2 <S> 4K S d<S> 



dcj) 2 



+ 



4 d 2 $ 



N 2 dz N 2 dz 



(7) 
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where N 2 = (K s — K) /F is the dimensionless Brunt- Vaisala frequency and we define $ = -d (Sp°/p) /dt. 

Equation (0 can be solved by separation of variables. The general solution that is regular as r — >• has the form 



&(r,(/),z,t) = F \ Jm(Ar) [A m cos(ra0) + £ m sin(ra</>)] Z m (z)T m (t) 



=o 



(8) 



where m ^ is an integer and A is determined by the boundary condit i ons. T he prefactor F is included as $ is expected to be 
of this order. This is the sa me result found by van Eysden &; Melatosl (|2008h but is slightly more general than the equivalent 
Abnev &; Eps tein as it allows for the possibility that FN 2 and K are of similar magnitude, a likely scenario in a 



neutron star (|van Eysden k, Melatos 2008). 



2.3 Boundary conditions 

The boundary conditions on $ are set by the boundary conditions on the O(E ) velocity fields, 

* - -mh (?) • 

as <I> is defined in terms of 5p° and is therefore O(E ) too. [To impose boundary conditions on the 0(E 1 ^ 2 ) flow, we would 
need to know dp 1 .] Assuming no penetration at the side wall, we have d^/dcj) = at r = 1 and hence 

oo oo 

<£(r, 0, Z,t) = F ^2 Jm(Xmnr) [Amn COs(?71(/>) + B rnn Sm(m(f>)] Zmn{z)Tmn{t) , (11) 

m—0 n—1 

where A mn is the nth root of J m (A) = 0. 

To find Z mn , we use the 0{E 1 ^ 2 ) axial flow, 

as Sv® — 0. We require 5v\ — at z — 0, so that the flow is symmetric about the midplane. The normalisation of Z mn is 
arbitrary, and we choose Z mn (l) = 1, giving 

(FN 2 - ft-) e?+ z - (FN 2 - P + ) e?- z 
Zmn{z) = (FN*-/3-)ef i +-(FN*-l3+)ef i - ' (13) 
with 

Another boundary condition applies to the top and bottom faces of the cylinder, which determines T mn . The mass flux 
into and out of the Ekman bounda r y layer at z = =bl is related to the circulation just outside this layer by (Pedlosky 1967; 
Walinl Il969l : lAbnev fe Epsteinl Il996l : Ivan Evsden fe MelatosJ l2008h 



K S ±(K 2 + N 2 \ 2 mn ) 1/2 ] . (14) 



*v*\ z =±i= t\e 1/2 [Vx{8v-v b )] z , (15) 

^ z=±l 

where is the dimensionless velocity of the boundary in the frame rotating at Q. Ekman pumping continues until the local 
fluid velocity, here Sv, matches the boundary velocity vb- For a rigid container, the final angular velocity equals Q + 5£1 in 
the inertial observer's frame, corresponding to = re^ in the rotating frame. To find T mn , we differentiate (|15p with respect 
to time and substitute equation JT2} into the left hand side of (p"5]) (note: 5v Q z = 0), and equations © and ([TO]) into the right 
hand side of JTSJ) . After some algebra, we find that the (ra,n)-th mode relaxes exponentially as $ oc exp(— cj mn t), with 

_ \ 2 mn [(FN 2 - p.) - (FN 2 - p + ) e"-] 



(4FK + \ 2 mn )(e^ -e?>-) 

Integrating $ with respect to time, the general solution for the pressure perturbation can be written as 

6p°(r,4>,z,t) 



(16) 



OO oo 



— C(r, 0, z) + F Jrn ( A mn r) ( A mn cos rncj) + B mn sin m<j)) Zmn (z)e " mnt , (17) 

m—0 n—1 

where A mn and B m n absorb a factor of u mn , and C(r,(j>,z) is the constant of integration. C(r,(/),z) is constrained by the 
boundary condition JTSJ) and must match the boundary velocity at z = 1. Using J9} and (fTU)) , we obtain 

"*rM,l) = - — ^ , (18) 
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VB + fafa l) = 



1 ac(r,<M) 

2P dr 



(19) 



2.4 Initial conditions 

All that remains is to specify the initial conditions, which determine A mn and B mn - Without specialising to a particular trigger 
for the spin- up event at t — or modelling the vortex unpinning and rearrangement that presumably accompanies it, we 
consider the general situation where these processes establish some instantaneously nonaxisymmetric p ressure field throughout 
the in terior. [Five possible physical causes of the nonaxisymmetry are discussed in detail in Section 1 of van Evsden &; Melatosl 
(2008).] We denote the initial state at t — by the symbol (*)P (r, 0, z) — Sp(r : 0, z, 0)/p(z). Specifying £P°(r, 0, z) is equivalent 
to specifying the initial velocity or density, which are related through J9]), (fTU)) , JT2}, and the O(E ) equation of motion, 



6p» 



d5p° 
dz 



(20) 



The choice of (*)P (r, 0, z) is arbitrary, but it should satisfy the boundary conditions outlined in Section [2.31 We eliminate 
C(r, 0, z) by evaluating (p"T[) at t = 0, obtaining 



oo oo 



= 8P°(r, tj>,z) + F ^ Jm(Xmnr) (A mn cosm<f> + B mn sinm^) Z mn (z) (e " m " t - l) 



m=0 n=l 

The coefficients Amn and Bmn are determined at 2 = 1 from £P°(r, 0, 2) and C(r, 0, 1) . In general, we have 
2 



A 



27T />1 

/ #/ drrJ m (A mn r)cos(m0) [&P°(r,<M) - C(r, 0,1)] 
^0 ^0 



(21) 



(22) 



P mn is given by the same formula, with cos(ra0) replaced by sin(ra</>). 



2.5 Velocity, density, and pressure solutions 

Equations ©, ([TO]) , (p"2]h ([20]) . and (|2T|) yield complete solutions for the velocity, density and pressure fields. Upon transforming 
back to dimensional variables and out of the rotating frame into the inertial observer's frame, we can write the results as 
follows: 



00 00 



(r,4>,z,t) = 5v r (r,(f),z,0) + -L 2 5fl ^ ^ — J m (A mn r/L) 



(FN 2 - f3-) e p + z/L - (FN 2 - £+) e 



P-z/L 



(FN 2 - P-) e^+ - (FN 2 - p+) e^- 
x {Amn sin[m(0 — fit)] — B mn cos[ra(0 — fit)]} (e~ E 1 UJrni%nt 



(23) 



V4, 



(r, <t>, z, t) = Qr + 5v+(r, </>, z, 0) + -LSQ A 

mn Jm(^mn^ / L) 



oo oo 



(FN 2 - 13-) e p + z ' L - (FN 2 - /?+) e p - z ' L 
(FN 2 - /?_) e"+ - (FN 2 - P+) e"- 



x {^ mn cos[m(</> - fit)] + B mn sin[m(c/> - fit)]} (e^ 1 ' 2 ""'""' - l) 



(24) 



OO oo 



V z (r,4>,Z,t) - -LSQE 1/2 ^ ^ \mnJrn(\mnr / ' L) 



d (3 + z/L _ e (3_z/L 

eP+ - eP- 



x {Amn cos[m(0 - fit)] + Bmn sin[ra(</> - fit)]} ^e E ' Wmnm - 1^ 



(25) 



p(r, 0, 0, t) = poe z ' Zs + <5p(r, 0, z, 0) + 



poLfldfl 



oo oo 



^ ^ ^ ^ Jm(^mnT / L) 



(FN 2 - p_) f3-e~P- z / L - (FN 2 - 0+) P+e'^ 1 



(FN 2 - P-) e^+ - (FN 2 - p+) e?- 



x {Amn cos[ra(0 - fit)] + Bmn sin[ra(0 - fit)]} ^e E 



(26) 



p(r, 0, t) = p gz s e z/Zs + 5p(r, 0, 0) + p L 2 fl5fl ^ ^ J m (A mn r/L) 



(PiV 2 - e-^- 2 ^ - (FN 2 - p+) e~P+ z/L 



(FN 2 - P-)eP+ - (FN 2 - p+) e^ 



x {A mn cos[m(0 - fit)] + Pmn sin[m(0 - fit)]} ^e s 7 UJrnnVtt - 1^ 



(27) 
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The initial velocity, density, and pressure are related to the chosen initial state £P°(r, 0, z), through, 

c / / p.\ 1 dSP°(r,<t>,z) f 
6v r {r,<l>,z,0) = - — (28) 

Svj, (r, )2 ;,O) = — — (29) 

d[p(z)5P°(r^,z)] 

6p(r, 0, z, 0) = — (30) 

Sp(r, 0, z, 0) = p(z)5P°(r, 0, 0) (31) 

In the limit t —> 00, an incompressible fluid spins up completely via Ekman pumping and approaches a steady-state 
solution, which matches the boundary at z = 1. In contrast, for a compressible, stratified fluid, part of the volume is 
untouched by Ekman pumping. In the latter case, the persistent, unaccelerated initial flow and the associated gradient in 
dissipate by viscous diffusion and adjust via inertial oscillations over the long time-scale 

2.6 5P°(r, 0, z) for a glitch 

In this paper, we assume that a glitch spins up the crust rigidly and axisymmetrically but that it initially excites nonaxisym- 
metric motions in the fluid interior; that is, 5v® and 5v^ are superpo sitions of cos(m0) and sin(m 0) modes immediately after 
the glitch. Possible physical mechanisms are outlined in Section 1 of van Eysden &; Melatosl (|2Q08h . The crust spins up rigidly 
to angular velocity Q + 5Q, which corresponds to C(r, 0, 1) = Fr 2 in equation (|22]h satisfying JT5J) and JT5J) as required. The 
arbitrary initial pressure perturbation 5P°(r, 0, z), which specifies the initial flow velocity through J9}, ([TO} and ([T2|) . is a sum 
of nonaxisymmetric modes satisfying the boundary conditions (e.g., no penetration of the side walls). In dimensionless form, 
in the rotating frame, we can write 
00 

6P°(r, 0, z) = F ^ C m r m (r 2 - l) cos(ra0) . (32) 

771=1 

No sin(ra0) terms or z dependence are included for simplicity, and the relative weights of the modes are parametrized by the 
constants Cm- We take C m = 1 for all m in this paper. 

The above initial condition is slightly more realistic than the one adopted by van Eysden &; Melatosl ( 2008h , who posited 



that the perturbed (spin-up) flow develops from 5v® — 5v® — i mmediately after the glit c h to a permanently nonaxisymmetric 
steady-state flow at the boundary [see equations (40) and (41) in Ivan Evsden fe Melatosl <|2QQ8h ]. There are two problems with 
the latter scenario. First, it involves nonaxisymmetric, and therefore nonrigid, motion of the top and bottom faces of the 
cylindrical container, which in reality would exert large stresses on the stellar crust, probably causing it to crack. Second, it 
artificially emits gravitational radiation in the steady state, even at t > E- 1 ^- 1 (cf. Section IO below). 

3 GRAVITATIONAL WAVE SIGNAL 

The gravitational radiation gen erated by the nonaxisymmetric spin-up flow in Section [2] is the sum of a mass quadrupole contri- 
bution, calculated previously bv lvan Evsden &; Melatosl (2008), and a current quadrupole contribution. The current quadrupole 
is typically smaller than the mass quadrupole by a factor ~ c/v. However, using the results of Section [2J the nonaxisymmetric 
velocity perturbation is larger than the density perturbation by a factor F, implying a wave-strain ratio /i mass /h CUTTe n t ~ 



Fc/v rsj Qc/g. We compute the current quadrupole wave strain in this paper and refer to Ivan Evsden &; Melatosl (|2008f ) for 
the mass quadrupole. 

3.1 Current quadrupole 

The far-field metric perturbation generated by a superposition of current multipole moments can be written as (Thorne 1980) 



00 1 



TT _ G V-^ V" d S m (t) B2,lm ;„i 

h » ~^dX. oti T f ' (33) 

1 = 2 m= — l 



in the transverse, traceless gauge, where t is the retarded time, D is the distance from source to observer, and T^ k 2,lm is 
a tensor spherical harmo nic which is a function of source orientation. The (Z,ra)-th multipole moment, S lm (t), is given by 



( Melatos fe Peraltallioiol ) 

1/2 



glm 327T 



(2Z + 1) 



: + 2 



2Z(Z — 1)(Z + 1) 



d 3 xr l x • cur\(pv)Y lm * , (34) 
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for a Newtonian source, where Y lm denotes the usual scalar spherical harmonic. In this paper, we only consider the leading 
order, quadrupole (I = 2) term. Importantly, S 2m depends only on the Fourier mode with frequency mQ in the spin-up 
flow described by equations (|23[) - ([27[) . In other words, the / = 2 metric perturbation is a linear superposition of terms each 
generated by a unique mode in the spin-up now. 

The plus and cross polarisations of the gravitational wave strain can be expressed compactly in terms of S 21 and S 22 . 
The axisymmetric Ekman flow leads to a quadrupole moment d 2 S 20 /dt 2 — O(E), which we neglect in this paper. Denoting 
the inclination angle between the rotation axis of the star and the observer's line of sight by z, we can write 

h+{t)= 2§D (^) 1/2 {MS 21 (t)}smi + Im[S 22 (t)}cosi} , (35) 

hx{t) = 4^D (^) 1/2 { Re [S 21 W]sin2i + Re[5 22 (t)](l + cos 2 i)} , (36) 
where an overdot symbolises differentiation with respect to time. 



3.2 Gravitational wave strain 



We compute the far- field metric perturbation at a hypothetical detector by combining the flow solutions in Section 12.51 with 
the boundary and initial conditions in Section 12.61 Appendix [A] shows how to rewrite the integral in (|34p to involve only 
the pressure perturbation 5p°, simplifying the evaluation of S lm . The final result for the current quadrupole moment, for 
< m ^ 2, takes the form 



n=l 



15m 



_ T/ \ -imQt , Tr -(E 1/2 uj rrin +im)nt 
mn v mn J c n - V mn^ 



with 

Umn — $ 



n,i / dr / < 
Jo Jo 



i t~t Tn ( 2 -i \ 

Ur ( r — 1 ) « 



V m n = dr dzr m+1 Z 2 ^UAmnJn 

Jo Jo 



(A mn r) 



(FN 2 - P-) e~P~ z - (FN 2 - p+) e _/3 + 2 
(FN 2 - (3-) e^+ - (FN 2 - £+) e?- 



where 

d 2 z d zm 2 d 2 . 
U = z— + r— — + 2F ( r 



dr 2 r dr 



drdz 



2 d 2 d 2 

— rz 

dz 2 drdz 



(37) 

(38) 
(39) 

(40) 



is a differential operator acting on everything to its right in equations ([38]) and ([39]) . Umn and V mn are straightforward to 
calculate analytically, but the full expressions are too lengthy to quote here. 

Substituting (|37p into (|35p and (|36p . we obtain the following expressions for the plus and cross polarisations as functions 
of time: 



h+(t) = h J2 



sin i j ({/in -Vin)smQt + Vi n e~ El/2uJlnnt [2E 1/2 u ln cos Qt - (Euo 2 ln - l) s'mQt] | 
-^cosi^4(U 2n -V 2 n)sm2Qt + V 2 ne~ El/2uJ2nnt [4E 1/2 uo 2n cos 2Qt - (Eu 2 2n -4) sin2Qt] 
sin2i{(Vi n - Uin) cosQt + Vi n e~ El/2uJlnnt [(Eu 2 n - l) cos Qt + 2E 1/2 u ln sin Qt] j 



(41) 



--(1 + cos" i){4(V 2n -U 2 n) cos 2Qt + V 2n e 



[(Eu 2n - 4) cos 2Qt + 4E 1/2 u 2n sin 2Qt] j 



with 



ho 



4ttGpqL 6 SQ Q 2 
3c 5 D 



(42) 



(43) 



Equations ([IT]) and ([42]) contain terms of order (Euo mn )° , (Euj^) 1 ^ 2 and (Eu^n) 1 • The derivation of the spin-up flow in 
Section [2] assumes E 1 ^ 2 <C 1. Over the range of values for K, N and E that we consider in Section [4] and [5] it is also true that 
E 1 / 2 ujmi <C 1. The quantity E 1 ^ 2 ujmn does become large for large n (uJmn — ► nit /2 as n —> oo), but the exponential suppresses 
the large- n terms and the infinite sum converges. For our purposes, truncating (|4ip and (|42)> at leading order O(E ) gives a 
good approximation. 
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In the scenario described in Section [2.61 the nonaxisymmetric initial perturbation is erased by Ekman pumping on the 
time-scale and the fluid spins up to rotate axisymmetrically with the boundary at z = 1. The effects of stratification 
and compressibility reduce the effectiveness of Ekman pumping, reducing the spin- up volume. As a result, some regions of 
the interior are incompletely spun up and preserve some of their initial nonaxisymmetric flow for t ^> unlike in the 
incompressible problem. The nonaxisy mmetry persists, emitting gravitational radiation c ontinuously, until viscous diffusion 
wipes it out on the time-scale ( Greenspan &; Howard! 19631 ; Benton &; Clark 19741 ). As the time-scale > 10 3 

years is comparable to, or greater than, the age of many glitching pulsars, one encounters the interesting possibility that 
neutron stars harbour a 'fossil' nonaxisymmetric flow in their interior, preserved by stratification, which continually emits 
gravitational radiation, and whose structure reflects the history of differential rotation and superfluid vortex rearrangement in 
the star. This possibility merits careful investigation in the future. It is not the same as the ar t ificial , nonaxisymmetric, nonrigid 
rotation of the crust postulated (for mathematical convenience) by Ivan Evsden fe Melatosl (|2008h (cf. also Section US). 



4 DETECTABILITY 

We now estimate the detectability of the gravitational wave signal derived in Section [3] by calculating the signal-to- noise ratio 
expected to be achieved by current- and next-generation long-baseline interferometers. The signal differs from a traditional 
continuous- wave source (e.g. an elliptical neutron star), because it decays over days to weeks (approximately 10 5 to 10 8 wave 
cycles). It is therefore counterproductive to integrate coherently past a certain time (if one ignores the fossil quadrupole 
discussed in Section |3?2|) . We find that the signal-to- noise ratio depends sensitively on the buoyancy, compressibility, and 
viscosity of the neutron star interior. For certain, plausible ranges of these variables, the signal is detectable in principle by 
Advanced LIGO. 



4.1 Signal-to- noise ratio 

The response of a laser interferometer to plus and cross polarisations h+(t) and h x (t) can be written as 

h(t) = F+(t)h+(t) + F x (t)h x (t) . (44) 



The b eam-pattern functions T + and F x depend on the rotation of the Earth and the sky position of the source (|Jaranowski et al.l 
19981 ). For the signal in Section [3l it is convenient to split h(t) into components that oscillate at the spin frequency of the 



star and its first harmonic, denoted h±(t) and h 2 (t) respectively. Writing h\, 2 (t) — i 7 +(t)/ii,2+(t) + Fx (t)hi^x (t) and keeping 
terms of order O(E ) in (gTJ and (g2]), we find 

oo 

h 1+ (t) = h sin i sin Qt fan ~ Vi n + Vine - ^ 2 " 1 " ') , (45) 

n=l 

oo 

h 2+ it) = -2h cos i sin 2ftt ^ (u 2n ~ V 2n + V^e^ 172 " 2 ^) , (46) 

71=1 
OO 

hi x (t) = -y sin2icosQ£^ (u ln - V ln + Vi n e~ El/ * ulnat> ) , (47) 

n=l 

oo 

h 2 x (t) = h (l + cos 2 i) cos 2Qt (^n - V 2n + y 2 ne" £l/2 ^ fit ) . (48) 

n=l 

The signal-to-noise ratio d for a quasi-dich romatic source (i.e . a sou rce consisting of two narrow-band peaks at frequencies 
/* and 2/*) is given by equations (80)-(82) in Jaranowski et al. ( 1998h . The result is 

T /2 2 /-To/2 

dt[/ii(t)] 2 + — — / dt[h 2 (t)] 2 . (49) 



Sh(f*) J-T /2 Sh(2f*j 7 -T /2 

In (|49|) . Sh(f) is the spectral noise density of the interferometer at frequency /, To denotes the total length of the coherent 
integration, and /* = Q/2tt is the stellar spin frequency. 

The integration time for a coherent search is normally limited by computational expense rather than the length of the data 
stream. Even when the radio ephemeris is known through radio observations, the radio and gravitational wav e phases may 



not b e equal, increasing the number of templates required for a search [e.g., the ^-statistic search for the Crab ([Abbott et al 



2008)]. We assume a computational limit of two weeks for the remainder of this paper. For the glitch recovery signal, the 
integration time is the minimum of the computational limit and the glitch recovery time-scale; integrating beyond the point 
where the signal decays away merely adds noise. The exact value of To which maximizes d depends on the search algorithm, 
but it is always of order the e _1 time constant for h(t), i.e. h(To)/h(0) = e _1 . For the general estimates below, we take 



© 2010 RAS, MNRAS 000.fTHl9l 



Continuous-wave gravitational radiation from pulsar glitch recovery 9 




decay time-scale of the leading (n — 1) term in equations (|45 [) -(|48 [) . The m = 2 mode decays 
1 mode, but the difference is moderate (1 ^ o^i/wii ^ 2) over the parameter space that we 



To = (£ 1/2 cj2ift)-\ the e 
more quickly than the m 
consider. 

Figure [T] illustrates how To depends on stellar parameters. The four panels in Figure [T] display contours of To (in days) 
on the K-N plane for four different values of E. The value of E in a neutron star is uncertain but Figure [T] demonstrates 



that it plays a significant role in determining To. One requires E 



1(T 



for the best match between (E ' cj2i^) and 



observed post-glitch recover y time- scales. This value is artificially lower than that expected from neutron- neutron scattering, 
E ~ 10" 7 (Q/rad s" 1 )" 1 ([Cutler k Lindbloml [l987l ; lAndersson et al.1 liooi : Ivan Evsden k Melatosl bold ) bec ause it is the 
effective value that arises when modelling the two-component Hall-Vinen-Bekarevich-Khalatn i kov superfluid (Per alt a et al. 
2005 ; Andersson k Com er 2006) as a single Newtonian fluid (Easson 19791 ; Abnev k Epstein 19961 ; Ivan Evsden k Melatos 



2008). 



To calculate d, we evaluate (|49|) with To = {E 1 ^ 2 uj2i^t)~ 1 and make several simplifying assumptions. First, we approximate 
h(t) by the leading (n — 1) terms in the infinite sums in (|45p ~ (|48p . For typical values of N and K, this introduces an error of 
< 10 per cent. Second, following I Jaranowski et all (| 19981 ), we average the functions sin(raOt) and cos(raQt), which oscillate 
much more rapidly than T + , F x , and exp(— t/To), over the observation period. The result is 

rT 



(1- 



S h (f* 



dt 



H — sin 2i t x 
4 



(1 



2 )h 2 A 2 (K,N) 



with 

Ai(K, N) 



1 



2e 



S h {2f.) 



V n )V a + ±Vl 



/•To 

/ <*[■ 
^0 



4cos 2 iT+ + (l + cos 2 i) 



(50) 



(51) 



l-e- 2V J l + e v 

As discussed in Section [3721 the signal is the sum of a persistent periodic signal associated with the fossil nonaxisymmetry 
(which decays on the long time-scale E^QT 1 ^> To) and the decaying signal generated by the Ekman flow. To be conservative, 
we only consider the latter signal, setting U m n = Vmn- Hence, (|51|) reduces to Ai(K,N) = 

The signal-to-noise ratio depends on the right ascension a, declination 5, and polarisation angle of the source as well as 
the location and orientation of the interferometer and the diurnal phase of the Earth. These quanti ties are usually known f or 



any specific source. However, to estimate detectability in general, we average d over a, 5, ip, and i ( Jaranowski et al. 1998) 
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Figure 2. Histogram of known glitching pulsars (dark shading) and observed glitches (light shading) as a function of frequency. Curves 
of anticipated spectral noise density for second- and third-generation interferometers are overlaid. [Note: glitches emit gravitational 
radiation at both /* and 2/* (see Section [3}.] Left panel displays Advanced LIGO configurations: zero detuning, high power (solid), black 
hole optimised (dashed), neutron star optimised (dash-dotted). Right panel displays ET configurations: conventional (solid), xylophone 
(dashed). 



i f 271 i r 1 i r 27r i r 1 



d(cosi) (. 



The beam pattern functions average to 



/ dtF* 
Jo 



/ AtFl 
Jo 



To . 2 > 



(52) 



(53) 



a, 5,ip WO / ck, 5,ip 

where £ is the angle between the arms of the detector. We give more details of this result in Appendix [Bj Substituting ([53]) 
into (|49|) and averaging over z, we obtain the following expression for the average signal-to- noise ratio: 

1/2 



Ai(K,N) + 4A 2 (K,N) 



Sh(f* 



(54) 



4.2 Second- and third-generation interferometers 



We now evaluate the signal-to-noise ratio (|54[) achieved by the second-generation interferometer LIGO, in both its Initial and 
Advanced configurations, and the third-generation, subterranean Einstein Telescope (ET). 

There are various detector configurations proposed for Advanced LIGC0- The best overall sensitivity across the entire 
frequency spectrum is achieved with zero detuning of the signal recycling mirror and high laser power. Below 40 Hz, the 
configuration optimised for 30M© black hole binary inspirals provides the best sensitivity. Above 40 Hz, the configuration 
optimised for 1.4M© neutron star binary inspirals provides the best sensitivity. However, the differences between the three 
configurations are small. 

Two configurations have been proposed for ET: a conventional interferometer ( Hild et al. 20081 ), and a dual-band xy- 
lophone config uration consis ting of two co-located interferometers, one optimised for low frequencies and the other for high 
frequencies (|Hild et alj|201(j ). Below 30 Hz, the xylophone configuration is more sensitive than the conventional configuration, 
by a factor of up to ~ 10 in the 5-10 Hz band. 

Figure [2] compares the spectral noise density of the different detector configurations. It also bins the number of known 
glitching pulsars and observed glitches as a function of frequency to illustrate which configurations are best suited for glitch 
searches. It is important to recall that the results of Section [3] predict gravitational radiation at both the pulsar frequency 
/* and 2/*, with the pulsar orientation determining which frequency has the stronger signal. The xylophone configuration 
of ET is the best choice for a glitch search. More glitches have been observed in objects with f < 30 Hz (83 per cent), 
where the xylophone is more sensitive, than / > 30 Hz (17 per cent) ( Per alt al 20061 : Melatos et al. 20081 ) . Additionally, the 
increase in sensitivity of the ET xylophone configuration over the conventional configuration is far greater below 30 Hz than 
the decrease above 30 Hz. In contrast, Advanced LIGO is not sensitive below 10 Hz and there is only a small difference in 
sensitivity between the different configurations over the frequency range where most glitches lie. As mentioned above, the 



1 LIGO Document Control Center: document number ligo-t0900288 -v3 
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K K 



Figure 3. Contours of angle- aver aged signal-to-noise ratio (d) versus normalised compressibility K and Brunt- Vaisala frequency N, for 
four existing and planned interferometers: (a) Initial LIGO, (b) Zero-detuning, high-power Advanced LIGO, (c) Neutron star optimised 
Advanced LIGO, (d) Black hole optimised Advanced LIGO, (e) Conventional ET, (f) Xylophone ET. Source parameters: /* = 100 
Hz, 5Q/Q = 2 x 10- 4 , E = 10~ 17 , D = 1 kpc. The detector spectral noise densities used are (S h (f*), S h (2f*)): Initial LIGO (1.75 x 
10- 45 Hz" 1 , 8.53 x 10- 46 Hz" 1 ), zero-detuning, high-power Advanced LIGO (1.59 x 10~ 47 Hz -1 , 1.39 x 10~ 47 Hz" 1 ), neutron-star- 
optimised Advanced LIGO (1.18 x 10~ 47 Hz -1 , 9.03 x 10~ 48 Hz -1 ), black-hole-optimised Advanced LIGO (3.77 x 10~ 47 Hz -1 , 1.84 x 
IO- 47 Hz" 1 ), conventional ET (6.68 x lO" 50 Hz" 1 , 6.68 x lO" 50 Hz" 1 ), and xylophone ET (1.56 x lO" 49 Hz -1 , 1.12 x lO" 49 Hz" 1 ). 



black-hole-optimised Advanced LIGO configuration is the most sensitive below 40 Hz but its advantage is slight and possibly 
outweighed by its slightly poorer performance at higher frequencies, where the strongest signals (from the fastest-spinning 
objects) arguably lie. 

Figure [3] displays contours of the average signal-to- noise ratio for Initial LIGO, Advanced LIGO (zero detuning and 
high laser power, neutron star optimised, and black hole optimised), and ET (conventional and xylophone) as a function of 
compressibility K and Brunt- Vaisala frequency N. The figure is produced for an object with /* = 100 Hz, E — 10" 17 , at a 
distance D — 1 kpc from Earth, with radius R = 10 km, mass M = 1.4M , po = 3M/4nR 3 , and g = GM/R 2 (Ekman pumping 
occurs in a thin surface layer, where g is uniform) . The step increase in angular velocity is taken to be 5Q/Q = 2 x 10 -4 , 
corresponding to the largest glitch observed to date ( Melatos et al]|2008h . The spectral noise densities [Sh(f*), Sh(2f*)] used for 
the six detector configurations are: Initial LIGO (1.74 x 10 -45 Hz -1 , 8.54 x 10 - 46 Hz -1 ), zero-detuning, high-power Advanced 
LIGO (1.59xl0 -47 Hz"\l.39xl0" 47 Hz" 1 ), neutron-star-optimised Advanced LIGO (1.18xl0 -47 Hz" 1 , 9.03 x 10" 48 Hz" 1 ), 
black-hole-optimised Advanced LIGO (3.77 x 10" 47 Hz" 1 , 1.84 x 10" 47 Hz" 1 ), conventional ET (6.68 x 10" 50 Hz" 1 , 6.68 x 
10" 50 Hz" 1 ), and xylophone ET (1.56 x 10" 49 Hz" 1 , 1.12 x 10" 49 Hz" 1 ). 

It is clear from Figure [3] that detectabil ity drops off sharply for K > 10. Buoyancy prevents Ekman pumping from 
spinning up the whole of the stellar interior ( van Evsden &; Melatosl 20081 ) . For large stratification (K s ^> FN 2 ), K s w K, 
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only a small volume of the interior is spun up and the current quadrupole is greatly reduced, with Ai,2 oc e~ 2K . There is little 
difference between the three Advanced LIGO configurations in panels (b), (c) and (d), or the two ET configurations displayed 
in panels (e) and (f) in Figure [3] All have similar sensitivity at 100 Hz. For ET, we find (d) > 3 for K < 10 and N < 1 and 
there is a reasonable possibility of detection. We require smaller values, e.g. N < 0.5 and K < 3, to achieve (d) > 3 with 
Advanced LIGO. 

To generalise the results in Figure [3] to an arbitrary object, we note that (d) scales with Ekman number as (d) oc E 1 ^ 4 
(square root of the number of cycles in the coherent integration). For relaxation time-scales of 3 to 300 days (Peralta 2006), 



and assuming K — N — 1, E ranges from 10 21 to 10 17 , which corresponds to an order of magnitude of variation in (d). 
The signal-to-noise ratio also scales with the spin parameters through the characteristic wave strain, viz. 



h Q = 6 x 10" 26 



( 10-") (l(/*Hz) (lkpc) • (55) 

The relative change in angular velocity SQ/Q is not necessarily equal to the observed glitch size 8v/v. In a vortex unpinning 
model, the two quantities are related through 5v/v ~ (I s / I c )(Ar / R)(5Q/£Y) , where I 3 /I c ~ 10 2 is the rat io of superfluid to 
crust moment of inertia , and Ar/R ~ 10 -6 is the normalised radial distance the unpinned vortices move (klpar et alj fl986: 



Melatos fe Peraltall2010h . Therefore, equating the observed glitch size to 5Q/Q yields a conservative estimate, given that SQ/Q 



may in fact be up to ~ 10 4 times larger. 

A coherent search synchronised to a radio ephemeris assumes that the radio and gravitational wave signals have the 
same phase. T his is not nece s sarily true. For example, in the landmark coherent J 7 - statistic search for the Crab pulsar in 
LIGO S5 data, [Abbott et "aL ( 20081 ) allowed for a fractional phase mismatch of up to 10 -4 . In our multiple scales analysis, 



we assume by construction that the nonaxisymmetric modes are stationary in the frame rotating with the pre-glitch angular 
velocity and remain so throughout the Ekman pumping process. In reality, the crust spins up to Q + and drags the 
axisymmetric part of the flow asymptotically to this increased angular velocity. Whether the angular velocity of the m ^ 
modes also increases during this process is unclear. It depends on exactly how the superfluid vortices repin following a glitch 
and rearrange themselves in a sheared Ekman flow, which is unknown at present. 

The number of templates required for a search can be estimated by modelling the frequency as /(£) = fo + fot. For 
a coherent search, the difference in phase between the model and gravitational-wave signals over the integration time must 
satisfy Acp < tt. For a two week integration, this corresponds to a maximum template spacing of Sfo = 3 x 10 -6 Hz and 
y =4x 10" 12 s~ 2 . 

During the glitch recovery, the frequency derivative is much larger than usual for an isolated pulsar spinning down 
electromagnetically. We approximate fo « Av/Tq. Conservatively, this yields fo ~ 10 -7 s -2 for a 100 Hz pulsar undergoing 
the largest glitch observed to date {Av/v ~ 10 -4 ) with an unusually short relaxation period of one day. This translates into 
a range of Afo — 10 -7 s -2 to search ove r and hence 3 x 10 4 templates in fo. To allow for some mismatch between the radio 
and gravitational wave phases we follow lAbbott et al. ( 20081 ) who searched over a window of ±6 x 10 -3 Hz centred on the 



radio frequency, i.e. Afo = 1.2 x 10 -2 Hz. Overall, therefore, a total of ~ 10 8 templates are required for a glitch search. 

The parameters N, K and E change the shape of the signal in two ways: the relaxation time is controlled predominantly 
by E, while the relative difference between the signals at /* and 2/* (in amplitude and relaxation time) is controlled by N 
and K. Our signal-to-noise ratio estimates in Figure [3] are based on the incoherent sum of the detector response at /* and 
2/*, so the relative phasing between /* and 2/* does not affect the detectability and the number of templates required. This 
would change in a more sophisticated search that combined the /* and 2/* responses coherently. 

To this point, we assume that radio observations provide the frequency, recovery time-scale, and trigger epoch for a glitch 
search. We now consider the scenario where this information is not known, as in a blind search. In the region of parameter 
space that we consider, 0.1 ^ N ^ 10, 0.1 ^ K ^ 10, and 10 -20 ^ E ^ 10 -8 , the minimum band width of the Fourier- 
transformed wave strain is w 6 x 10 -12 /*. Hence, searching over the frequency range 1-600 Hz requires ~ 10 12 templates in 
fo multiplied by ~ 10 4 templates in fo as discussed above. 

In addition, the sky positio n, time of occurren ce, and recovery time are unknown for a blind search. In a LIGO search 



for unknown periodic sources ( Abbott et al. 2003), the sky is divided into 31500 patches. The lack of an electromagnetic 



trigger means that the data must be searched in many blocks, starting, for example, one day apart (coherent integration over 
a shorter recovery time is unlikely to be detectable) and integrating over increasing lengths of time, up to the computational 
limit, to account for the fact that To is unknown. A proper estimate of the associated computational expense lies outside the 
scope of this paper. 



5 CONSTITUTIVE PROPERTIES OF BULK NUCLEAR MATTER 

Figure [3] clearly demonstrates that the strength of the gravitational wave signal depends sensitively on the constitutive 
properties of bulk nuclear matter (e.g., the equation of state) and its dissipative or transport coefficients (e.g., viscosity). We 
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Figure 4. Ratio of fundamental and first-harmonic Fourier amplitudes |/i_|_(/*)|/|/i_|_(2/*)| (dashed blue contours) and full widths at 
half maximum r_|_(/*)/r_|_(2/*) (solid red contours) for six slices of the four-dimensional parameter space (N, K, E,i). In each panel, 
two variables are fixed, with K = 1, N = 1, E = 10 -17 , and i = 7r/4 as appropriate. 



show that these properties can be inferred in principle from the detailed shape of the gravitational wave signal. The results of 
this approach can be linked to terrestrial experiments, e.g. with heavy ion colliders, although there is an important distinction 
between ~ GeV collisions of ~ 10 2 nucleons in a terrestrial particle accelerator and ~ 10 57 static nucleons at ~ MeV energies 
in a neutron star. 

In a real search, one seeks to extract param eters like K and N by fitting a template to the interferometer data in the time 
domain (Clark et al . 2007; Havama et al.l l2008). However, to illustrate the scientific potential of the fitting exercise, we Fourier 
transform h+(t) and h x (t) and focus on the gross features of the spectrum. We neglect the permanent fossil quadrupole (see 
Section |3?2|) and assume that there is no interference between peaks. The four peak amplitudes and four peak widths (at /* 
and 2/*) of h+(f) and h x (f) provide enough information to solve for K, N, E, i, and ho by matching to the theoretical 
predictions in (|45p ~ (|48p . We take ratios to eliminate ho (which depends on the unknowns po, R, SQ, and D) and focus on the 
remaining parameters. 

Figure [4] displays six slices through the four-dimensional parameter space. Contours are shown for the amplitude ratio 
|/i+(/*)|/|/i+(2/*)| and the width ratio r+(/*)/r+(2/*), where r+ jX (/*) is the full width at half maximum of the peak in 
\h+ iX (f)\ centred at /*. The figure is drawn for the parameter ranges 0.1 < K ^ 10, 0.1 ^ N ^ 10, 10~ 20 < E < 10~ 8 , and 
^ i ^ 7r. We evaluate the first 20 terms in the infinite sums in (|45p ~ (|48p . In those panels where K, N, E, and i are held 
fixed, we use the fiducial values K — 1, N = 1, E — 10 -17 , and % — ty/4 respectively. 

The inclination angle determines the relative strength of the m — 1 and m — 2 modes of /i+ and h x through the tensor 
spherical harmonic T^ k 2,lm in (|33p . The contours of |/i+(/*)|/|/i+(2/*)| are nearly vertical in panels (d), (e) and (f) of Figured 
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In fact, if we consider additional amplitude ratios, we can infer i independently from the other parameters. Dividing the Fourier 
transforms of (g5) by (|i7|), and (46]) by (|I8|), we obtain \h+(f*)\/\h x (/*)! = seci and \h+(2f*)\/\h x (2f*)\ = 2 cosi/(l + cos 2 i) 
respectively. These expressions overdetermine z, yielding its value and an independent cross-check. The inclination angle can 
also be inferred from the radio or gamma-ray pulse profile and polarisation swing by assuming a par ticular emission model 
( Lvne h Manchester! [l988l : Iffibschman fc AronslbOQll : feai h Spitkovskvll2009UChung h Melatosll2010l ). The width ratios are 
independent of i\ the time-scale over which the signal decays does not depend on the location of the observer. This is illustrated 
in panels (d), (e) and (f) of Figure H] where the contours of r+(/*)/r+(2/*) are horizontal. 

The compressibility K and Brunt- Vaisala frequency N are inextricably linked in the sense that they feed into both the 
amplitude and width ratios in a complicated manner. However, once we determine i according to the formula above, we can 
immediately extract K and N from panel (a) of Figure U as the value of E does not influence any of the amplitude or width 
ratios (see below). By plotting contours of the measured ratios of |/i+(/*)|/|/i+(2/*)| and r+(/*)/r+(2/*) on the K-N plane, 
we can read off the values of N and K from the intersection point of the contours. One might be tempted to use other 
amplitude and width ratios as a cross-check on K and N, or to break the degeneracy in the case of multiple intersection 
points. However, most of the ratios are related trivially through the inclination angle and supply no additional information, 
e.g. r x (/.) = r + (/.),r x (2/,) = r+(2/,), |fc x (/+)| = cosi|M/+)l, and |ft x (2/+)| = (cosi + seci)|M2/ + )|/2 

The Ekman number E is important in determining the recovery time-scale and hence the Fourier width. It also appears 
in ho through To. However, it influences all peaks in the same way and drops out of all amplitude and width ratios. In panels 
(b), (c) and (f) of Figure U] the amplitude and width ratio contours are vertical. As mentioned in Section [4] an approximate 
value of E can be inferred from the e- folding time of h(t), as K and N only weakly influence this quantity. However, if K 
and N are known, e.g. by following the procedure described in the above paragraph, we can determine E from the absolute 
peak widths. Finally, ho can be determined from the absolute peak amplitudes once the values of all the other parameters are 
known. 

Future gravitational-wave measurements of the compressibility, viscosity and Brunt- Vaisala frequency of bulk nuclear 
matter can be compared to a range of terrestrial experiments and theoretical calculations. The compressibility is commonly 
expressed in terms of the compression modulus k, which is relate d to our normalised compr essibility through K = Am p gR/K,, 
where A is the mean atomic numbe r and m v the proton mass (van Eysden &; M elatos 2008). He avy-ion collisions and nuclear 



resonance experiments measure k ( Sturm et al. 2001 ; Vretenar et al.1 20031 ; Piekarewiczl 20041 ; Hartnack et al. 20061 ) . Com 
pressibility can also be infer red from the symmetry energy measured in hea vy- ion collisions or obtained through neutron-skin 
thickness measurements ( Chen et al. 1 120051 : iLi et al.l I2OO8I : IXu et al.1 [2009) . The shear viscosity is often expressed in terms 
of the ratio rj/s, where s is the specific entropy. It is related to the Ekman number by E = (A 'ks I 'm v R 2 Q)(ri / ' s) where 
1 ^ A' ^ 2 is the entropy per nucleon in uni ts of Boltzmann's constant (Ivan Eysden fc Melatos! 20081 ) . The shear viscosity has 
also been measured in heavy- ion collisions ( Adler et al. 2003 ; Adare et al.l 120071 ). Neutron stars are stably stratified because 
the c oncentration of charged particles increases with density but chemical equilibrium is maintained (Reisenegg er &; Goldreich 
1992). Stratifica tion provides a buoyancy force proportional to the Brunt- V aisala frequency squared, which has been calculated 



theoretically ( Reisenegger &; Goldreich 1992 ; Lai 1994 ; Passamonti et al. 20091 ) 

In Table [T] we quote a selection of experimental and theoretical values for K, N, and E under neutron star conditions. 
Dimensionless values of N and E assume 0/2tt = 100 Hz. In line 1, th e compression modulus is inferred from the ratio of 
the K + multiplicity in Au+Au and C+C collisions at ~ GeV energies ( Sturm et al. 2001 ; Hartnack et al. 20061 ) . In line 2, 
the compression modulus is obtained b y fitting a relativist ic mean- field mode l to the distribution of isoscalar monopole and 
isovector dipole strengths of Zr and Pb (Vreten ar et al . 2003; |Piekarewiczll2004l ). In line 3, the compressio n modulus is obtaine d 
from the measured nuclear symmetry energy from isospin diffusion in heavy- ion collisions ( Chen et al. 20051 ; Li et al. 2008) 



Line 4 lists the ratio of s hear viscosity to specific entropy measured in A u+Au collisions at an energy of 200 GeV (Adle r et al 



20031 ; Adare et al. 2007 ). Theoretical calculations of shear viscosity by ICutler fc Lindbloml (jl987l ) for neutron-neutron and 



electron-electron scattering, corresponding to the normal and superfluid states respectively, are listed in lines 5 and 6. More 
exotic states, which ma y exist in the neutron star core, will have a different viscosity. Line 7 lists the shear viscosity due to 
quark-quark scattering ( Jaikumar et al. 20081 ) . In lines 8 and 9, we quote calculated values for t he Brunt-Vaisala frequency, 
the l atter including centrifugal forces in a rapidly rotating star ( Reisenegger &; Goldreich! 19921 ; Lai! 19941 ; Passamonti et al. 
20091) . 



6 CONCLUSIONS 

In this paper, we calculate analytically the gravitational radiation emitted during the post-glitch recovery phase by the 
nonaxisymmetric Ekman flow excited by a glitch. The calculation is done in the context of an idealised, cylindrical star with 
a uniform viscosity, compressibility, and stratification length-scale. We compute the signal-to-noise ratio for current- and 
next-generation long-baseline interferometers and find the following promising result: for a large glitch (50/0 = 10 -4 ) from 
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Table 1. Experimental and theoretical results for compressibility, viscosity and Brunt- Vaisala frequency. 



Quantity 


Experiment /Theory (E/T) 




Result 


Dimensionless 


Reference 


K 


Au+Au and C+C collisions (~ GeV) (E) 




k w 200 MeV 


K = 0.97 


1,2 




nuclear resonances (E) 




k « 240-270 MeV 


K = 0.72-0.81 


3,4 




nuclear symmetry energy (E) 




k = 210 MeV 


K = 0.93 


5,6 


E 


Au+Au collisions (200 GeV) (E) 




f]/s « U/AnkB 


E = 8x 10- 20 


7,8 




neutron-neutron scattering (T) 


V 


= 2 x 10 20 g cm" 1 s _1 


E = 5 x 10- 9 


9 




electron-electron scattering (T) 


V 


= 6 x 10 20 g cm" 1 s _1 


E = 1 x 10~ 8 


9 




quark-quark scattering (T) 


V 


= 5 x 10 15 g cm -1 s _1 


E = lx 10- 13 


10 


N 


chemical composition (T) 




AT* ~ 500 s- 1 


N = 0.8 


11,12 




centrifugal correction (T) 




N = 0.32-0.84 


N = 0.32-0.84 


13 



(I) Sturm et al. (2001) , (2) Hartnack et al. (2006) . (3) Vret enar et al. (2003) . (4) Piekarewicz (2004), (5) Chen et al. (2005) . 
(6) Li et al. (2008) . (7) Adler et al. (2003) . (8) Adare et al. (2007) . (9) Cutler fc Lindblom (1987) . (10) Jaikumar et al. (2008), 

(II) Reisenegger fe Goldreich (1992) , (12) Lai (1994) , (13) Passamonti et al. (2009) 



a neutron star D — 1 kpc from Earth and spinning at /* = 100 Hz, the angle- averaged signal-to-noise ratio (d) exceeds three 
for N < 0.5, K < 10, and E ~ 10" 17 with Advanced LIGO and N < 1, K < 10, and E ~ 10" 17 with ET. 

Perhaps the most obvious shortcoming of our idealised model is its cylindrical geometry. There is a n oble history of 
using a cylinder to model spherical astronomical objects and also in classical geophysical studies of the Earth ([Pedloskvll 19671 ; 



Wali n 1969; A bnev fe Epsteinll 19961 : Ivan Evsden fe Melatosl 2008), because it admits analytic solutions, which in general have 



not 



yet been found for a sphere. We ignore magnetic fields for simplicity, although they are large in neutron stars ((Cutler 
2002|), interact with the superfluid ( Mendelll 1998 ), and therefore modify Ekman pumping. We model the interior of a neutron 



star as a single Navier-Stokes fluid, whereas in reality it is a multi-component superfluid, consisting o f superfluid neutrons 
and superconducting protons wh ich interact with each other via mutual friction and entrainment (e.g. Lattimer &; Prakashl 
20041 : lAndersson &; Comer! 120061 ) . The spin- up process in a coupled multi-component fluid of this kind, in the presence of 



gravitational stratification and compressibility, is an unsolved and difficult problem. 

In our model, the crust accelerates instantaneously from Q to Q + SQ and remains at this higher angular velocity. A 
more realistic model would co nserve total angular momentu m by solving self-consistently for the response of the crust to the 
viscous back-reaction torque ( van Evsden &; Melatos 201p| ). In the context of the present model, we can approximate this 
effect crudely by replacing the glitch size SQ/Q at t = with the permanent frequency jump after the recovery ceases. None 
of the conclusions change qualitatively. 

Understanding the glitc h mechanism remai ns an unsolved problem. Glitch waiting times are exponentially distributed and 
their sizes fit a power law (Melat os et al . 2008), indicative of inhomogeneous collective behaviour on large scales, e.g. vortex 
avalanches. In contrast, nuclear structure calculations suggest that the area densit y of pinning sites (e.g. lattice defects) is much 
greater than the area density of vortices (|jonesll2002l : Donati &; Pizzocherol 2003 ), suggesting that the system is homogeneous 
on large scales (pinned Abrikosov array). The gravitational wave signal calculated here helps to discriminate between these 
two views, as it is a measure of the internal nonaxisymmetry. From a simple, random walk argument, the largest relative 



glitch size that arises from vortex movement in a star containing n vortices is SQ/Q 



-1/2 



. If the value of SQ/Q inferred 



from a gravitational- wave detection approaches this maximum, it is safe to infer that large-scale inhomogeneities are present. 
Note that we take C m = 1 in Section T2.6I However, if only a fraction of the internal flow is nonaxisymmetric, C m should be 
reduced in proportion. In vortex unpinning models SQ/Q can be up to four orders of magnitude larger than the observed 
glitch size (see Section H~2)) , leaving considerable scope to get detectable gravitational- wave signals. 

Vortex unpinning theories of glitches rely on the build up of a lag between the crust, which spins down electromagnetically, 
and the superfluid, whose rotation is fixed by the number of vortices, until a glitch is triggered. We know that the lag does not 
disappear completely after the glitch (i.e. co-rota tion is not restored) because a reservoir effect (i.e. glitch size oc waiting time) 
is not observed in glitch data ()Wong et al.ll200ll ); only a small, random fraction of the lag relaxes during a single event, and 
that fraction is determined by the microscopic history of the system, as in any avalanche process. In this model, we assume 
conservatively that the crust and fluid co-rotate before the glitch. However, in the more realistic scenario just described, there 
is ongoing differential rotation between crust and core, suggesting that glitching pulsars may continuously emit gravitational 
radiation. 

Another possibility leading to a continuous gravitational wave signal beyond just the post-glitch recovery period is the 
'fossil flow' discussed in Section [3] Stratification prevents Ekman pumping from spinning up the whole interior, leaving a 
remnant of the initial nonaxisymmetric flow untouched. This flow emits gravitational radiation until damped over the much 
longer diffusion time-scale. If so, we may be able to extend the coherent integration time beyond the recovery time-scale, 
increasing the likelihood of detection. Even more intriguing is the possibility that any neutron star which has experienced 
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differential rotation in its past retains some part of this fossil flow for > 10 3 years, thereby bearing an imprint of the star's 
formation and rotation history. We plan to study the matter fully in a following paper. 

For a typical neutron star at a distance of 1 kpc, the signal-to-noise calculations in Section U argue that there is a 
reasonable chance interferometers like Advanced LIGO or ET will detect the largest glitches. The outlook is more optimistic 
if we consider nearby 'dark' neutron stars. For the estimated galactic population of ~ 10 9 neutron stars (cf., ~ 1800 radio 
pulsa rs discovered to date), recent Monte-Carlo simulations predict the closest objects are located ~ 8 pc from Earth (jOfekl 
2009). At this distance, Initial LIGO is able to detect the largest glitches with (d) > 3 for N < 1 and K < 10 and Advanced 
LIGO is sensitive to smaller glitches with SQ/Q > 10 -6 . However, the signal frequency, glitch epoch, and sky position are 
unknown electromagnet ically, so searching for 'dark' glitches is a difficult proposition. None the less, our results suggest 
cautious optimism about the chances of detecting a glitching (or otherwise differentially rotating) neutron star with the next 
generation of gravitational- wave interferometers. 
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APPENDIX A: SIMPLIFYING x • CURL(pv) 

It is straightforward to evaluate S lrn by substituting (|23p ~ (|27p directly into (|34|) . However, the calculation is easier and 
more transparent if we first simplify the integrand in f|34|) to depend only on 5p. Expanding according to p \— > p° + 5p, 
v H> v° + Sv = rVLecf) + Sv, we express the integrand to first order as 

x • V X (p°v° + 5pv° + p°Sv) . (Al) 

The first term in (|A1[) is independent of time. It does not emit gravitational radiation, so we discard it. The second term in 
(|A1[) reads 

d 



V X (5pv°) = Q [iz + rz4~ - r 2 -^-) Sp , 
\ or dzj 



dz 2 



d 2 



drdz 



(A2) 
(A3) 



To move from JA2J) to (|A3| we use the Navier- Stokes equation JTJ to first order in Rossby number and zeroth order in Ekman 
number. In the rotating frame and neglecting the centrifugal term, as in Section [2] it reads 

2p°(ft x Sv) = -V5p - g5p , 

from which we obtain 



1 d6p 

x~ ■ 

9 dz 

The third term in (|A1[) can be rewritten in a similar way. From (|A4[h we find 
2Qp°e z X (e z X Sv) = V X (Spe z ) , 

where (e r , e<p, e z ) are the basis vectors in cylindrical coordinates. Noting that 5v° z — 0, as the axial flow is 0(E 1 ^ 2 ) 
left with 

P°Sv = -^V X (6pe z ) , 
and the third term in (|A1[) is 



(A4) 
(A5) 



x • V X (p°<Sv) 



2Q 



drdz 



'dz* 



Sp 



Combining (|A3[) and (|A8|) , and replacing d 2 /d(j) 2 by — m 2 , we arrive at 



x • V X (pv) 



2Q 



d 2 z d 

z 1 

dr 2 r dr 



zm 



drdz 



dz 2 



drdz 



• 2z- 



dz 



6p° 



(AG) 



(A7) 



(A8) 



(A9) 



There is a subtle issue around neglecting the centrifugal correction to (|A4[) , which is of order F. If we evaluate S lm by 
substituting (|23]) - ([27)) directly into (|3l|) . we implicitly include centrifugal terms in x • V X (pv) (by virtue of failing to exclude 
them explicitly). This approach is internally inconsistent, because centrifugal terms of this order are excluded from the flow 
fields (|23[) - ([27[) following the assumption in Section f2.2l leading to J5} and ©. It is therefore preferable to evaluate (|34|) for 
S lm using (|A9[) , so that the centrifugal correction to the zeroth-order structure is consistently excluded from both the flow 
fields and S lm . 



APPENDIX B: BEAM PATTERN FUNCTIONS 

The complete expressions for the beam pattern functions are ( Jaranowski et al. 1998), 
F+ (£) = sin C[a(t) cos 2ip + b(t) sin 2ip] , 
F x (t) = sin C[b(t) cos 2ip - a(t) sin 2ip] , 
with 

a(t) = ^ sin 27(3 — cos 2A) (3 — cos 25) cos[2(a — cj) r — Vt r t)\ — i cos 27 sin A(3 — cos 25) sin[2(a — <j) r — £l r t)] 

1 13 

+ — sin 27 sin 2 A sin 25 cos [a — <j) r — £l r t] — — cos 27 cos A sin 25 sin [a — <j) r — Q r t] + — sin 27 cos A cos 5 

b(t) = cos 27 sin A sin 5 cos[2(a — (f) r — Q r t)] + i sin 27(3 — cos 2A) sin 5 sin[2(a — <j) r — Vt r t)\ 



+ cos 27 cos A cos 5 cos [a — <j) r — £l r t] + — sin 27 sin 2 A cos 5 sin [a 



(Bl) 
(B2) 



(B3) 



■ n r i] . (B4) 
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The right ascension and declination of the gravitational wave source are given by a and 5 respectively, and ip is the polarisation 
angle. The latitude of the detector is denoted by A, Q r is the angular velocity of the Earth, and <j) r is the diurnal phase of the 
Earth. The angle counterclockwise between East and the bisector of the i nterferometer arms is 7 , and the angle between the 
arms of the interferometer is £. We average over a, 5 and ip according to ( Jaranowski et al. 19981 ) 

{■■■) aM = ^zl daxi/ d(sin<5) x i- / (...) • (B5) 



We evaluate (J^° dtF+) a ,s,ip and (J* dtF 2 ) a ,$^ for use in Section [H Averaging (|B1|) and (|B2[) over ip, we obtain 



T 



dtFl 



T 



dtFi 



1 . 2, 
-sm C 



To 



dt([a(t)] 2 + [b(t)} 2 



(B6) 



All the dependence on a and S is contained in a(t) and b(t). After some straightforward but lengthy algebra, we find that the 
dependence on all other angles drops out, leaving 

<[a(t)] 2 + [6(t)] a ) a , = |. 

Substituting (|B7[) into ()B6|) and evaluating the now trivial time integration we obtain the result stated in equation (|53[L 



(B7) 



rTo 



dtFl 



To 



To . 2 > 
-sm C- 



(B8) 
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